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A novel Thomas-Fermi (TF) approach to inhomogeneous superfluid 
Fermi-systems is presented and shown that it works well also in cases 
where the Local Density Approximation (LDA) breaks down. The nov- 
elty lies in the fact that the semiclassical approximation is applied to 
the pairing matrix elements not implying a local version of the chemical 
potential as with LDA. Applications will be given to the generic fact 
that if a fermionic superfluid in the BCS regime overflows from a nar- 
row container into a much wider one, pairing is substantially reduced at 
the overflow point. Two examples pertinent to the physics of the outer 
crust of neutron stars and superfluid fermionic atoms in traps will be 
presented. The TF results will be compared to quantal and LDA ones. 



1. Introduction 

The quantal treatment of pairing in inhomogeneous systems is a notoriously 
difficult problem. This is especially true for systems containing a large 
number N of particles as it is usually the case for cold atoms in traps 
(N ~ 10 6 |3or even for smaller systems if they are deformed as can happen 
for nuclei. Semiclassical approaches may be very helpful in such cases. The 
simple and very well known Local Density Approximation (LDA]P is not 
always applicable because for its validity the condition that the size of the 
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Cooper pair (coherence length), £, must be smaller than a typical length 
I over which the mean field potential is varying (I is, e.g., the oscillator 
length in the case of a harmonic potential) is not always fullfillcd. We here, 
therefore, will apply the TF approximation directly to the pairing matrix 
elements whose evaluation only requires the usual TF condition that the 
wave lengths involved must be smaller than . We will show that, indeed, 
our approach also works for cases where £ is larger than I where the LDA 
fails. We will demonstrate this for the BCS approach in this paper. 

The physical systems we are interested in concern cold atoms in traps 
and nuclei, both in so called overflow or drip configurations. For the lat- 
ter overflow or drip means that there is such a large neutron excess that 
the selfconsistent mean field container is full up to the edge. In the inner 
crust of neutron stars where the nuclei form a Coulomb crystal these extra 
neutrons overflow into the interstitial space and form there a more or less 
dense neutron gas which also can be superfluid. In the inner crust the nu- 
clei can actually turn into sheets and the neutron gas can form in between 
the sheets (a so-called lasagne configuration®). As a first example we will 
treat in a schematic model such a slab configuration as is shown in Fig. rjj 
mostly because the quantal solution is readily available and, therefore, can 
serve as a test case for the validity of the TF approximation for treating the 
pairing problem. Indeed, we will find that the TF approach reproduces the 
quantal solution of the pairing properties, besides some shell fluctuations, 
very accurately. On the physical side, we point out that at the overflow 
point pairing can be strongly suppressed. This finding will then also be 
reproduced with a system where cold atoms are filled into a spherical con- 
tainer consisting of a narrow part at low filling, suddenly going over into a 
much wider container at higher chemical potentials, as it is displayed in the 
left panel of Fig. [3f . A slightly different situation occurs with a double 
well potential, as the one shown in the left panel of Fig. [5] . Again this 
potential is used in a slab configuration and TF and quantal results for the 
gaps are compared. 

At the end, we return to nuclei in the crust of neutron stars where they 
are embedded in a more or less dense gas of neutrons. A Wigner-Seitz cell 
approximation will be applied to investigate this situation. Again similar 
features as in the previous examples will be found around the overflow point 
in the transition from the outer to the inner crust. 
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Fig. 1. Left: Schematic view of the slab used in this work with a perspective view of 
the potential which is translationally invariant in x, y direction. Right: Quantal and TF 
pairing gaps in the slab geometry as a function of the chemical potential. 

2. The formalism 

In this section, we will explain the TF approach to the pairing problem and 
apply it first to the slab configuration shown in Fig. Q] Let us start out by 
writing down the usual BCS equations in three dimensions 



where the e^s are the single particle energies and fi the chemical potential 
which can be used to fix the particle number N — with 



The wave functions and eigenenergies of a box as shown in Fig. Q] with 
a potential-hole are given irP . For pairing, we use a contact force with 
a cut off A, to make things simple. The single particle states in a slab 
configuration then become \v) — |n,p > where n are the discrete quantum 
numbers in transverse direction and p the momentum quantum numbers 
in slab direction. To obtain the gap equation in this case we start by 
integrating the gap equation (fj]) over momenta in slab direction: 




(1) 
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with E n (p) = yj (e n + s p — p) 2 + A 2 the quasiparticle energy, 9(x) the step 
function, and £ n ,e p being the discrete single particle energies in transverse 
direction and kinetic energies in slab direction, respectively. After simple 
algebra, one arrives at the following gap equation for a slab configuration 

A ™ = - Yj Q ( A - £ n>)V nn >K n ,. (4) 
n' 

The pairing tensor in equation (4) is then given by 



*r - m a i A-/x+ y/(A-^) 2 + A2 

K » - ln £n -, + v^-^+a^ (5) 

where m is the particle mass and the indices n stand for the level quantum 
numbers in the confining potential of the left panel of Fig. [1] The matrix 
elements V nn ' = —g \(p n (z)\ 2 \cp n >(z)\ 2 dz of the pairing contact force 
v P air(r — r') = —gS(r — r') used in this case can be evaluated straightfor- 
wardly from the wave functions (p n [z) given irP . 

Before we show the results, let us explain our Thomas-Fermi (TF) ap- 
proach for this problem. In the weak coupling regime, we have A//i << 1. 
In this case the canonical basis-^ can be replaced by the Hartree-Fock or 
mean-field one: 



H\n) = e n \n). (6) 

At equilibrium and for time reversal invariant systems canonical conjuga- 
tion and time reversal operation are related by 



(r|ra) = (n\r) =$> (r^nn) = (ri|p„|r 2 ), (7) 
with p n = \n)(n\. For the pairing matrix element, we, therefore, can write 

Vnn' = {nn\v\n'n') = J (r 2 |/5„|ri)(rir2|w|rir' 2 )(r'i| / 5„|r' 2 )dri(ir2dr' 1 (ir' 2 . 

(8) 

The Schroedinger equation (j6|) can be writen in terms of p n as 

(H - e n )pn = 0. (9) 

Taking the Wigner transform of this latter equation, we obtain in the h — > 
limit the following c-number equatiorf : (H c i, — e)/ c (R, p) = 0. The 



March 2, 2013 17:43 World Scientific Review Volume - 9in x 6in 50yearsl0 



Using World Scientific 's Review Volume Document Style 5 



solution of this equation in the sense of distribution theory is with x8{x) = 
given by 



MR,p) 



9 TF (E) 



S(E-H cL ) + 0(h 2 ). 



(10) 



with 



H A = ^L m+V (R) and 9 ™(E) = j dKd^E - H cl ). 

with m*(R) the effective mass and V(R) the mean field potential. Equation 
(JTUJ) means that the phase space distribution corresponding to a state \n) at 
high energy is concentrated around the classical energy shell that, indeed, 
is a well known fact. 

The TF version of the gap equation (jj) then reads 



,A 

A(E) = - dE'g(E')V{E,E')K(E') (11) 

JVo 

with K(E) an obvious generalisation of K n in ([5]). The matrix elements 
V(E,E') can be evaluated in replacing \ip n (z)\ 2 by 3 



= / - WW) S^>" 2|E - v " r 1,2 ■ (12) 

which is the on-shell TF density in transverse direction (please note that 
we are in a ID case here, contrary to what is treated above where it is 3D). 
As the reader will easily realise, the way of proceeding is very different from 
usual LDA where the finite size dependence is put into the (local) chemical 
potential, /i(z) = [/, — V(z), whereas here it is put into the matrix elements 
of the pairing force (notice that in LDA they are computed using plane 
wave functions). 



3. Results 



We are now in a position to solve the quantal and TF gap equations in 
the slab geometry for the confining potential displayed in the left panel of 
Fig. [TJ As an example we take as cut off A = 50 MeV counted from the 
edge of the pocket potential whose depth be Vb = —40 MeV. Its extension 
ranges from —R to +R with R = 10 fm. The wide potential with infinitly 
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Fig. 2. (Coloronline) Left panel: Position dependence of the gap in the slab geometry 
for different values of the chemical potential. Quantal, TF, and LDA results are shown. 
Notice that A for /i = 0.5 and -5.0 MeV is practically zero in the gas region. Right 
panel: Comparison of quantal (dots) and TF (broken lines) values of the pairing tensor 
K. 




- co /(o -0.2 

_ mag opt 


U/lico o t = 


84.34 


"CO /(I) =0.15 

mag opt 


A/h(0 o i a 


164.34" 


ra mllJ /(o opt - 0. 1 f 'Z- s \ 






■ f X 














1 'if 

1 \ '.y 

l\ // 





20 40 60 80 100 120 140 160 181 
uVhcfl 



Fig. 3. (Coloronline) Average TF gaps at the Fermi energy as a function of the chemical 
potential (right panel) for the potential shown in the left panel. In the completely filled 
optical trap (fj, = U) we accomodate 10 5 atoms in each spin state. The total number of 
atoms in the trap with fi/hw pt =40, 80 and 120 are indicated in the upper horizontal 
axis. 



high walls has extension from — L to +L with L = 100 fm. For the coupling 
strength we take g= 150 MeV fm 3 . The result for the gap at the chemical 
potential /i is shown in the right panel of Fig. [T] as a function of /i. We 
start with /z from the bottom of the pocket well, i.e. with zero density. 
We then increase /i, i.e. the density. We see that once the fill up of the 
pocket reaches its top, the values of the gap sharply drop and practically 
reach zero. In the continuum the gaps slowly rise again. We see that 
quantal and TF values are in close agreement. The overshoot of the TF 
solution for negative /i very likely is due to the smallness of the pocket 
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Fig. 4. (Coloronlino) Pairing energy as a function of the chemical potential for the slab 
(left panel) and the double harmonic potential discussed in the text (right panel). 





Fig. 5. (Coloronline) Double well potential in slab geometry (left) and gaps as a function 
of the chemical potential (right). The different curves correspond to the parameters a=6 
and b=4, a=8 and b=4 and a=8 and b=2 (see text for explanation). Please note that 
TF and quantal values cannot be distiguished on the scale of the figure, besides in the 
region of the dips where TF passes through the average of the quantal oscillations. 

which only can accomodate nine bound levels. The bunches of resonances 
in the continuum of the quantal solution are interesting but we did not try 
to explain them in this work. Before we come to an explanation of the 
drop of the gaps at overflow (drip), let us study the gaps as a function of 
position in transverse direction. Quantally the position dependent gap is 
defined as: A(z) = —gK(z) with K(z) = "^2, K n \Lp n (z)\ 2 . Semiclassically, 
the relation between the gap and the pairing tensor becomes: 



A(z) = -g f dEg TF (E)K{E)p T /(z 



(13) 
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Fig. 6. (Coloronline) Left: Binding energies per particle as a function of average density 
in a Wigner-Seitz cell. Red dots indicate quantal Skyrme HF calculations by Negele- 
Vautherin and black dots correspond to semiclassical results with the variational Wigner- 
Kirkwood (VWK) method and the Gogny D1S force. Right: BCP equation of state for 
the inner crust compared with the predictions of the Baym-Bethe-Pethick one. 
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Fig. 7. (Coloronline) Left panel: Radial dependence of the TF gaps in the considered 
WS cells. The end points indicate the radius of the WS cells. Right panel: Gap values 
for the Wigner-Seitz cells corresponding to Fig. [6] The zero point indicates the transition 
from isolated nuclei (black dots) to the dripped case (red diamonds). 



In the left panel of Fig. [21 we show the gap profiles for three different 
values of the chemical potential: /j, = 40, 0.5, and - 5 MeV. We see that 
quantal and TF results agree, up to shell fluctuations, very well. We also 
show the LDA results. They can be locally as wrong as by 50 percent. For 
other choices of system parameters the LDA error may even be worse. This 
stems from the fact that in TF (and, of course, also quantally), there is 
coupling between inside and outside the pocket, i.e. the Cooper pair wave 
function extends into both regions what tends to equilibrate the values of 
the gaps. In LDA the contrast is much too strong. The drop of the gaps 
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Fig. 8. (Coloronline) Comparison between TF and LDA gaps as a function of the po- 
sition in a WS cell containing a single |g°Zr nucleus (left) and gg 00 Sr (right). 
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Fig. 9. (Coloronline) Left panel: The TF (blue diamonds) and LDA (green triangles) 
gaps at the origin are compared with the gaps of the free neutron gas (A(iJ ce ;;), red 
dots) as a function of the Fermi momentum of the neutron gas at the edge of the cell. 
Right panel: pairing energy per nucleon in the inner crust. The TF results are compared 
with the quantal HFB values of ref.HH 



when crossing the threshold can be explained by the fact that the single 
particle states are strongly delocalised in the outer container and, thus, 
their contribution to the pairing matrix element V n ^ n i becomes very small. 
In the right panel of Fig. [5] we show the quantal and TF pairing tensors, 
K n and K(E) respectively, defined before. We emphasize again the close 
agreement between quantal and TF results. 

Having gained faith into our TF approach, we now can explore other 
geometries and other systems, which are more difficult for quantal solu- 
tions. In the right panel of Fig. [3] we display the result for the gap A in 
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Fig. 10. (Coloronlinc) Neutron pairing gaps averaged with the pairing tensor (uv) along 
the Z=40 and Z=70 isotopic chains obtained with the Shlomo potential as mean field 
and the finite range Gogny force renormalized by an attentuation factor of 0.85 in the 
pairing channel. Filled dots correspond to the BCS calculation and solid thick lines to 
the TF approach. 





Fig. 11. (Coloronline) Left: TF (dashed line) and fluctuating (solid line) level den- 
sity (upper panel) and accumulated level density (lower panel). See text for details. 
Right: Averaged gaps along the Sn isotopic chain computed fully quantally at BCS level 
using the Gogny D1S force (circles), semiclassically using the fluctuating level density 
(diamonds), and with the TF approach (solid line). 



the spherical double harmonic oscillator potential shown in the left panel of 
Fig. [3] The latter may be realised with cold fermionic atoms to study the 
overflow situation. A zero range pairing force with strength g—-1.0 and cut 
off A=164.34 (in the corresponding optical trap units with oj opt = 2n x 1000 
Hz taken froirP) is used in this case. We see that the result is qualitatively 
similar to the slab case, though in this spherical geometry the dip does not 
quite reach zero and also is shifted slightly to an energy above the break 
point of the potential. Note that this depends strongly on the choice of the 
ratio tOmag/uopt as it can be seen in the figure. Also the gap starts to de- 
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crease towards the minimum quite early. This is contrary to what happens 
in the slab case, where the change is very abrupt. The reason probably 
lying in the spherical symmetry of the considered system. It would be in- 
teresting to see whether our prediction can be verified experimentally. 

It also is interesting to study the pairing energy, defined in TF ap- 
proximation as E pa i r — j? J dEg(E)A(E)n(E). The results of the pairing 
energy per particle for the slab potential, Fig. [TJ and in the H.O. case of 
Fig. [3] are shown in Fig. 2J We see that E pa i r behaves quite differently in 
the two cases. In the spherical example for cold atoms the depression at 
the overflow point is also seen in the pairing energy whereas the depression 
is completely washed out in the case of the slab. The reason for this quali- 
tative and strong difference must come from the fact that in the slab case 
the drop of the gap as a function of the chemical potential A(/j,) at overflow 
is extremely steep, almost vertical. Furthermore, in the pairing energy cor- 
responding to the slab, the pairing tensor k(E) should actually be replaced 
by K(E) corresponding to <j5j> - Being integrated over the momenta in slab 
direction it does not show any peak at E = \l as is the case in the spherical 
case. Therefore, in the integral of E pa i r also gap values further away from 
the overflow point are picked up which are not small at all. 

Another interesting geometry which can be considered is a potential 
with a barrier at the origin, i.e. a double well potential, see the left hand 
panel of Fig. [5] We treat this exemple again in slab geometry as in Fig. 
[TJ It roughly may mock up an oblate and very elongated trap potential 
for cold atoms with a double well potential in transverse direction. The 
results for the gaps A(//) are shown in the right hand panel of Fig. [5] This 
potential, which depends on two parameters a and b, is defined as follows: 
V{z) = \{z + a) 2 i£z<b, V(z) = + ±(1 - f )z 2 when \z\ < b and 

V(z) = |(z — a) 2 for z > b. As before, TF and quantal results are in 
excellent agreement. When the chemical potential reaches the top of the 
barrier, there occurs again a reduction of the gaps, since the wave func- 
tions at the Fermi surface suddenly get more extended above the barrier. 
This is the same effect as in the previous examples, although it is in this 
case less pronounced. Such a double well potential allows for the creation 
of a Josephson current if the population in left and right well are out of 
balance. 9 Our TF approach may strongly facilitate the description of this 
phenomenon in the case of cold fermionic atoms. 
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Let us now make a more realistic study of Wigner-Seitz (WS) cells to 
simulate the inner crust of neutron stars. In this approach one considers 
a single nucleus of N neutrons and Z protons inside a spherical box of 
radius R ce u as well as a uniform background of Z electrons to preserve 
the charge neutrality of the celP . The mean- field, as explained iiP^ ( 
is computed selfconsistently in the TF approach using the BCP energy 
density functional^ . In this semiclassical calculation we consider the same 
WS cells and mass numbers as in the old quantal calculation of Negele and 
Vautherirf^ . However, as far as shell corrections are not included, in our 
semiclassical calculation, we take as representative nucleus in each cell, the 
beta-stable one computed a la TF along the corresponding isotopic chain. 
This is why the atomic numbers Z of the representative nuclei differ from 
the ones reported irP^ while their mass numbers A coincide. 

It must be pointed out that the total energy per baryon obtained with 
our TF approach is in very good agreement with the quantal values reported 
irP , as it is explicitly discussed in Reffl^ and again shown here in the left 
panel of Fig. [6l As an additional test of our TF mean-field calculation, in 
the right panel of Fig. [6l we also display the EOS (i.e. pressure as a function 
of the WS average density) in the inner crust obtained in our semiclassical 
calculation compared with the results provided by the Baym-Bethe-Pethick 
EOS^ which is considered a benchmark in large scale neutron star calcu- 
lations. We find excellent agreement between both calculations. 

The semiclassical description of the WS cells including pairing corre- 
lations at TF level is obtained from this mean-field using the finite range 
part of the Gogny D1S forced in the pairing channel . In the left panel 
of Fig. [Jj we display the radial dependences of the gaps in some selected 
WS cells. It is seen that when the gap is small outside the region of the 
nucleus, then the gap also is small inside the nucleus. This stems from the 
very large coherence length where one neutron of a Cooper pair can be in 
the huge volume of the gas and the other inside the small volume of the 
nucleus (proximity effect). In this way the gas imprints its behavior for the 
gap also inside the nucleus. Such a conclusion was also given in a quan- 
tal Hartree-Fock-Bogoliubov (HFB) calculation by Grasso et al. irP^ what 
shows that the here employed BCS approximation apparently yields very 
similar answers as a full HFB calculation for WS cells^"^ . More precisely, 
let us point out that the gaps in the region of the nuclei, corresponding to 
the inner crust and displayed in this figure, are strongly affected by the 
neutron gas. To illustrate this fact, we display in the right panel of Fig. [7J 
the values of A(R = 0) (blue diamonds) and A(i? = R ce ii) (red filled cir- 
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cles), compared with the gaps of the free neutron gas (continuous black 
line) at the density corresponding to the edge of the cell p(R ce u). The 
semiclassical TF gaps A(R = R ce ii), as expected, closely follow the free 
neutron gas values in agreement with HFB calculations^ . As seen, the 
gap values at the origin, A(R — 0), are also strongly correlated with the 
gaps of the free neutron gas. For small average densities below p ~ 0.02 
fm~ 3 , the A(R = 0) values are larger than the corresponding gaps at the 
edge of the WS cell, as it can also be appreciated in the left panel of the 
figure. However, for larger values of the average density in the WS cell, this 
tendency is reversed and the gap at the edge is larger than at the origin 
pointing to the increasing influence of the neutron gas. These conclusions 
can also be drawn from an independent quantal study in RefPS where in 
the right panel of Fig. 4 very similar features for the local gap values in 
different WS cells can be seen as in our TF study. We also show in Fig. 
7, right panel, the LDA values of A(R = 0) by the filled (green) triangles. 
We remark that those values undershoot quite strongly the corresponding 
TF values at the higher densities. 

For further illustration of this effect, we show in the two panels of Fig. [8] 
a comparison between LDA and present TF results for the gaps in two 
particular WS cells. In the case of the largest cell whose representative 
nucleus is 4Q°Zr, we see locally a huge difference in the surface region of 
the nucleus. This simply stems from the fact that in this case the gap 
is very small and, therefore, the coherence length very large invalidating 
LDA. A study with examples a little less unfavorable for LDA is given irP^l 
. This wrong behaviour of the local LDA gaps at low average densities 
can also be seen in Fig. 6 of RefP^ . From that figure we conclude 
that our TF calculation reproduces, qualitatively, the global trends of the 
quantal gaps at low average densities. The behaviour of the semiclassical 
gaps at high average density is clearly different and it is dominated by 
the neutron gas as it can be seen in the right panel of Fig. |8] where the 
gap of the representative nucleus || 00 Sr is plotted as a function of the 
radius. Locally, LDA and TF show a depression in the center and the 
gap increases with increasing distance till it reaches its neutron gas value. 
The central depression is stronger by about 30 percent in LDA than in 
TF. This behaviour is similar to the one exhibited by the quantal gaps 
compared with the LDA ones displayed in RefP^ . In other words, this 
means that the contrast between inside and outside of the nucleus is much 
too pronounced in LDA, however, quantally as well as in TF this contrast 
is strongly attenuated by the proximity effect. The semiclassical TF gaps 
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at the Fermi level are displayed in the left panel of Fig. [H] as a function 
of chemical potential \i. In the inner crust, i.e. for positive values of 
H, they show a similar behaviour as the gaps at the edge of the WS cell 
displayed in the right panel of Fig. [7J This behaviour can be expected 
as far as the gap at the Fermi level is rather an average quantity and, 
therefore, strongly influenced by the neutron gas as it also happens in the 
HFB calculations of Ref ™ . In this figure we also include the gaps of some 
WS cells corresponding to the outer crust where all neutrons are bound 
with negative values of the chemical potential /i. Again we see that the 
gap practically vanishes at zero chemical potential when the neutrons start 
to drip. In the right panel of Fig. |H] we display the pairing energies per 
nucleon corresponding to the WS cell of the inner crust. These energies 
are also correlated with the neutron gas and display a similar behaviour 
as the one of the gaps at the Fermi energy. Again the pairing energy per 
nucleon vanishes when neutrons arrive at the drip. It is rewarding that in 
the low average density regime the TF pairing energies per nucleon follow 
the same trend as the quantal HFB ones 21 shown by the (red) filled triangles 
(however, a slightly different model with a somewhat stronger pairing force 
than in our case is used there). 

For isolated nuclei at the neutron drip the situation may be somewhat 
different. First, it may be that in this situation the difference between 
HFB and BCS approaches is more significant. Also strong shell fluctuations 
surely play an important role. Somewhat conflicting results in this respect 
exist in the literature. In reP^ very similar results to ours are found for 
iS-wave pairing, see Fig. 4 of this reference and also the discussion about 
it irP^ . On the other hand irPH g a p seems to rise towards the drip 
before it bends down. Similar results have also been found irJ22J . The HFB 
calculation of Hamamoto has recently been repeated and extended passing 
from negative chemical potentials to positive ones and it was found that 
the S'-wave gap clearly continues down to zero, touching zero at a slightly 
positive value of /lP^I _ j n explaining why in other works the gap is rising 
towards the drip, one has to keep in mind that an average gap should 
be calculated with the pairing tensor and not with the density matrix. 
The latter picks up the gaps at all energies which may not be small at 
all, even though the gap at the Fermi level is very small, see the right 
panel of Fig. [T] On the other hand an average with the pairing tensor 
generally only picks up the (small) gaps around the Fermi level. It also 
is intuitively clear that for other than S'-waves gaps the situation will be 
somewhat different. This is due to the finite centrifugal barrier which keeps 
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the wave function concentrated on the domain of the nucleus as long as the 
corresponding energy stays below the barrier. However, large scale HFB 
calculations around the neutron drips of nuclei indicate that in general 
pairing is reduced at the drip lineP^ In Fig. 10 we show a schematic study 
which may qualitatively reflect the real situation. A Z, N dependent Woods 
Saxon potential (without spin orbit) given by ShlomcE^ was taken as the 
mean field and the BCS equation has been solved with the Gogny D1S 
pairing forc d 14 * 1 ^ . Isotopic chains for two values of Z have been claculated. 
For Z = 70 the drip practically coincides with a shell closure of the neutrons 
and, therefore, the gap falls to zero at the drip for this case. On the other 
hand, for Z = 40 the neutron drip does not coincide with a closed shell and 
then the gap has substantial values around the drip. Globally, however, 
a clear decreasing tendency of the gap towards the neutron drip can be 
observed as also reflected by the TF values. The fact that neutron gaps 
decrease with increasing isospin was actually pointed out long time ago, 
see Refs.L^U^ _ Real nuclei at the neutron drip may be either spherical 
or deformed (see RefP^). For spherical drip nuclei it often happens that 
neutrons are at or very close to shell closure whereas for deformed nuclei 
this is not the case. The two situations then resemble the scenario displayed 
in the right and left panels of Fig. [10] though, there, we imposed sphericity 
also for the case Z = 40. 

A very promising possibility to recover shell effects is given by the fact 
that the latter mostly stem from the shell effects in the level density. We 
have discussed this problem with some detail in an earlier publication^ . 
The basic idea is to replace in the TF gap equation ((6|) the semiclassical 
g{E) by its quantal counterpart, slightly smeared out by using gaussians 
centered at the quantal eigenvalues so that one obtains a continuous func- 
tion of E without destroying substantially the shell structure as is shown in 
Fig. [TT] Inserting this into the gap equation (j6]), allows to recover almost 
completely the full quantal gaps as is shown in the right panel of Fig. [TT] 

In this work, we exclusively have treated pairing in BCS approximation. 
However, for certain situations and quantities, the more general HFB ap- 
proach may be mandatory as, e.g., in the cases of rotation or a magnetic 
field. It is relatively straightforward to generalise the BCS-TF approach 
also to the HFB case. For this one has to consider fully non diagonal ma- 
trix elements {nin2\v\n\in2') ■ In the matrix elements, we have replaced 
\ip n (z)\ 2 by the TF expresion for the on shell density p^ F (z). In the off 
diagonal pairing matrix element, we need wave functions and not densities. 
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Therefore, in TF approximation, we then can use <p n (z) — * yP~E F ( z ) f° r 
the individual wave functions. Of course, this complicates the solution of 
the gap equation but this is always the price to pay when passing from BCS 
to HFB, quantally as well as semiclassically. The TF-HFB approach shall 
be investigated in future worlP^ . Let us finally mention that for spherical 
systems the TF approach can be generalised to partial waves as was done 
for the pairing matrix elements irP. 

4. Summary 

Summarizing, we have studied superfluid fermions in a large container, ei- 
ther external (cold atoms) or created self consistently (nuclei) for situations 
where the top of the fluid reaches the edge of a small pocket located at the 
origin of the wide confining potential. The gap drops to zero at the edge 
before rising again when the density fills up the outer container. This at 
first somewhat surprising phenomenon can be explained quite straightfor- 
wardly. Such situations, as already mentioned, can exist in cold atoms and 
nuclei in the inner crust of neutron stars, two examples treated here with 
their specific form of containers. For small systems, like isolated nuclei at 
the neutron drip, the situation may be blurred by shell effects. 

As an important second aspect of this work, we showed that a novel 
Thomas-Fermi approach to inhomogeneous situations can cope with situ- 
ations where LDA fails. This means that our TF approach is free of the 
restrictive condition, prevailing for LDA, that the Cooper pair coherence 
length must be shorter than a typical length I (the oscillator length in the 
case of a harmonic container) over which the mean field varies appreciably. 
On the contrary, our TF theory has the usual TF validity criterion, namely 
that local wavelengths must be shorter than I. 

The accuracy of our TF approach opens wide perspectives for a treat- 
ment of inhomogeneous superfluid Fermi-systems with a great number of 
particles not accessible for a quantal solution of the BCS (HFB) equations. 
Such systems may be cold atoms in deformed containers (eventually reach- 
ing millions of particles), superfluid- normal fluid (SN) interfaces, vortex 
profiles, etc. As a matter of fact, as is well knowrP , the TF approach 
becomes the more accurate, the larger the system. Thus the TF approxi- 
mation is complementary to the quantal one in the sense that the former 
works where the latter is difficult or even impossible to be obtained numer- 
ically. 
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Ideas and part of this paper have been published in earlier works, see 
for instance RefsP^^ . A similar semiclassical approach also has been put 
forward for mesoscopic systems in RefP^ . 
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to publication. Special thanks are due to A. Pastore and J. Margueron 
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supported by the IN2P3-CAICYT collaboration (ACI-10-000592). One of 
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